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Cancer cells utilize large amounts of ATP to sustain growth, relying primarily on non-oxidative, fermentative 
pathways for its production. In many types of cancers this leads, even in the presence of oxygen, to the secretion 
of carbon equivalents (usually in the form of lactate) in the cell’s surroundings, a feature known as the Warburg 
effect. While the molecular basis of this phenomenon are still to be elucidated, it is clear that the spilling 
of energy resources contributes to creating a peculiar microenvironment for tumors, possibly characterized by 
a degree of toxicity. This suggests that mechanisms for recycling the fermentation products (e.g. a lactate 
shuttle) may be active, effectively inducing a mutually beneficial metabolic coupling between aberrant and non- 
aberrant cells. Here we analyze this scenario through a large-scale in silico metabolic model of interacting 
human cells. By going beyond the cell-autonomous description, we show that elementary physico-chemical 
constraints indeed favor the establishment of such a coupling under very broad conditions. The characterization 
we obtained by tuning the aberrant cell’s demand for ATP, amino-acids and fatty acids and/or the imbalance 
in nutrient partitioning provides quantitative support to the idea that synergistic multi-cell effects play a central 
role in cancer sustainment. 


INTRODUCTION 

At heart, a cell’s energetic problem consists in selecting 
how to process nutrients (say, glucose molecules) into chem¬ 
ical energy (adenosine 5’-triphosphate, ATP) that will then 
be transduced into useful forms of mechanical or chemical 
work. Rapid cellular growth, in specific, requires high rates 
of macromolecular biosynthesis and of energy production, 
which presupposes (a) fast ATP generation, and (b) tight con¬ 
trol of the cell’s redox state, i.e. that the ratio between the 
levels of electron donors and acceptors stays in a range that 
guarantees functionality. Most often, molecular oxygen is the 
primary electron acceptor in cells, playing a central role in the 
electron transfer chain (ETC) that constitutes the main ATP- 
producing mechanism in cells. When a glucose molecule en¬ 
ters the cell, it is normally metabolized by glycolysis, a highly 
conserved reaction pathway that converts each glucose anaer¬ 
obically into two molecules of pyruvate, with the concomi¬ 
tant production of 2 ATPs. In presence of oxygen, cells can 
operate the ETC, which begins with the conversion of pyru¬ 
vate into acetyl-coenzyme-A (acetyl-CoA). The reaction path¬ 
ways responsible for the subsequent production of ATP (and 
of many macromolecular precursors like amino-acids) are the 
Tricarboxylic Acid (TCA) cycle and Oxidative Phosphoryla¬ 
tion (OXPHOS). These complex groups of reactions (roughly 
100 processes altogether in the bacterium E. coli) are able to 
generate the largest energy yield in terms of molecules of ATP 
produced per glucose molecule intaken (up to 36, adding to 
the 2 given by glycolysis), and release carbon dioxide as a 
waste product. In absence of oxygen, however, cells cannot 
rely on the ETC and the ATP yield of glycolysis (2) is to a 
good approximation all the energy they can generate. In such 
conditions, the pyruvate obtained from glycolysis is then re- 
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duced to other carbon compounds (e.g. acetate, ethanol, lac¬ 
tate) that are normally excreted in variable amounts. The con¬ 
version of pyruvate to lactate, is carried out by a single reac¬ 
tion catalyzed by the enzyme lactate dehydrogenase (LDH). 

The energy-generating strategies just described are, in a 
sense, the two extremes, and cells usually operate mixtures 
of the two even in presence of oxygen, leading to ATP yields 
below the theoretical maximum of 38 (typically around 30). 
However, fast proliferating cells normally display high rates 
of glucose intake and produce ATP anaerobically even in the 
presence of oxygen, thereby spilling potentially useful carbon 
and energy resources. A hint about why a large glucose in¬ 
flux may favor the use of lower-yield pathways is provided 
by the fact that processing high glucose fluxes via glycolysis 
requires high rates of production of adenosine 5’-diphosphate 
(ADP) and of NAD+, via oxidation of NADH. The simplest 
way to convert NADH back into NAD+ is by reduction of 
pyruvate to lactate via LDH. Therefore sustaining high rates 
of glucose metabolization may imply lactate overflow. This 
however seems to suggest that a cell with a large glucose in¬ 
take should always prefer to generate energy by glycolysis. 
Therefore, different constraints (physical, regulatory, thermo¬ 
dynamic, etc.) may be at work in the selection of a cell’s ener¬ 
getic strategy Cl. We note that recent high-throughput studies 
of the compounds secreted by growing bacteria in controlled 
environments (the so-called exo-metabolome) uncovered that, 
besides the standard outputs of overflow metabolism, a pre¬ 
viously unsuspected diversity of molecules accompanies the 
excretion of carbon equivalents El- 

Likewise, aerobic glycolysis with lactate overflow (a.k.a. 
Warburg effect) is found to occur in many types of cancers 
EH, although it cannot be considered as a necessary signa¬ 
ture of malignancy 0. In order to explain its predominance 
at the molecular level, several ideas have been pursued, from 
structural or genetic abnormalities in the mitochondria 0, to 
the roles played by the hypoxia-inducible factor (HIE) 13 and 
by a specific isoform of the glycolytic enzyme pyruvate ki- 
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nase il, to seemingly unrelated genetic events leading to the 
increase of glucose transporters and the loss of growth con¬ 
trol in cancer cells M- The search for a molecular basis of 
the Warburg effect has given momentum to the idea that treat¬ 
ments targeting the simpler energy-producing apparatus of a 
cell (as opposed to its complex and not entirely known ge¬ 
netic profile) carry a higher potential than gene-based thera¬ 
pies CQl. Unluckily, our understanding of cancer metabolism 
is still largely incomplete iniini, and research on the role 
of metabolism in tumor genesis and progression is currently 
undergoing a major revival na. 

Cell-autonomous models of cellular metabolism based on 
genome-scale reconstructions of the underlying network of 
metabolic reactions have been widely studied, proving effec¬ 
tive in explaining the origin of overflow metabolism in differ¬ 
ent contexts. Constraint-based approaches like Flux-Balance 
Analysis (FBA), for instance, have shed light on the roles that 
different types of constraints may play in driving the selection 
of energetic strategies towards aerobic glycolysis. An interest¬ 
ing and intuitively appealing suggestion is put forward by the 
macromolecular crowding scenario m, according to which 
enzymes carrying out energy-efficient pathways can at most 
occupy a fraction of the intracellular volume (e.g. the mito¬ 
chondria in eukaryotes). Because reactions are assumed to be 
enzyme-limited, the activity of aerobic pathways will neces¬ 
sarily level out once all deputed enzymes work at full speed, 
while that of glycolysis may still increase. These models are 
indeed able to predict carbon excretion by unicellulars ca 
as well as the Warburg effect in cancer cells lfT4l [Thl . Many 
issues however deserve further analysis. 

One in particular concerns the fact that cellular waste prod¬ 
ucts such as lactate tend to deteriorate the extracellular envi¬ 
ronment, so that the sustainability of aerobic glycolysis de¬ 
pends crucially on the possibility that such a pollution is re¬ 
mediated within cancer’s microenvironment. This raises ques¬ 
tions about the role of tumor-stroma metabolic interactions. 
Tumors are universally characterized by a marked upregula- 
tion of glucose transporters mini, which allows them to 
largely surpass their non-aberrant (stromal) neighbors in the 
competition for nutrients. As a consequence, stromal cells 
must adjust their energetics in ways that have only started to 
be analyzed d). 

Lactate exchanges are increasingly being recognized to 
play an important role in this respect. In a scenario that 
has been characterized by genetic signatures over the past 
few years, tumor-derived lactate can be taken up by glucose- 
starved cancer-associated fibroblasts (CAFs) m, which 
would then both help their own survival and foster cancer’s 
aberrant growth by sanitizing the environment. Here we 
present a quantitative in silica analysis of the above mecha¬ 
nism. Our results confirm that tumor-to-stroma lactate shut¬ 
tling appears as a robust consequence of basic physical and 
chemical constraints, in a manner that is largely independent 
of (i) the metabolic function (e.g. energy production, amino- 
acid synthesis, fatty acid synthesis, etc.) that aberrant cells 
strive to maximize and, perhaps more surprisingly, (ii) of the 
degree to which they manage to maximize it. The core of 
our conclusions can already be reached through a highly sim¬ 


plified model of metabolism that captures the bare essential 
features of the mechanism. The intuition is then validated on 
a large-scale model of catabolism based on recent reconstruc¬ 
tions of the human reactome. After characterizing the network 
at the cell-autonomous level, we consider a system formed by 
two such networks sharing the same glucose supply, focusing 
on the case in which one cell aims to maximize a growth- 
related objective function (e.g. energy production) with the 
other subject to a simple ATP maintenance constraint. The 
emerging scenario is studied by tuning the ‘degree of max¬ 
imization’ by the aberrant cell through a rigorous sampling 
scheme that provides a statistically significant description of 
all feasible metabolic phenotypes compatible with the con¬ 
straints. Lactate overflow and functional tumor-stroma cou¬ 
pling appears robustly upon maximizing production of ATP 
as well as of several biomass precursors. 

It is important to stress that a tumor-to-stroma coupling, 
while compatible with the observation that the vast majority 
of cancer cells display glycolytic metabolism, is not the only 
lactate-based shuttling postulated to be relevant for tumor pro¬ 
gression (see e.g. 1201). Other mechanisms and their relation 
to the current study are outlined in the Discussion. 

The remainder of the manuscript is organized as follows. 
The Results section begins with a brief summary of a styl¬ 
ized single-cell model of aerobic glycolysis presented in ifT^ . 
where it is shown how a crowding constraint can divert ATP 
production from oxidative phosphorylation to fermentation. 
The simple model is then extended to two cells whereby the 
extent of lactate shuttle is computed analytically. A large scale 
metabolic network of the human core metabolism (HCCN) is 
then introduced. A rigorous sampling of the metabolisms pro¬ 
duced by HCCN confirms the possibility of cell coupling by 
means of lactate shuttle between a cancer cell and a stromal 
cell. The rigorous samplings allows us to compute the corre¬ 
lation coefficients among the fluxes of the cancer and stromal 
cell. In the Discussion section, experimental evidences for 
lactate shuttle are presented and. Anally, physiological impli¬ 
cations and possible experimental validation of the results are 
discussed. 


RESULTS 


Metabolic fluxes and crowding constraint in a minimal 
cell-antonomons model 

Our starting point is the highly simplified cell-autonomous 
model defined in El, in which lactate secretion is related 
to an ad hoc constraint limiting volume available for energy 
production in cells that maximize the ATP output flux. In par¬ 
ticular, when ATP demands exceed the volume allocated for 
ATP production, the cell produces ATP also through fermen¬ 
tation. In the model of Vazquez et al., a cell is characterized 
by a glucose intake flux a glycolytic flux /giyc (yield¬ 

ing two pyruvate molecules per intaken glucose) and an over- 
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all mitochondrial OXPHOS flux /ox that transforms pyruvate 
into energy (a table summarizing the most relevant variables 
defined here is given as Table SI). Alternately, pyruvate can 
be turned into lactate by lactate dehydrogenase (LDH), with 
flux /ldh- Lactate can then be expelled from the cell so as 
to avoid acidification. Mass balance and steady state condi¬ 
tions impose that metabolite levels do not change in time, so 
that, for example, the net production and consumption fluxes 
of pyruvate must compensate, implying 

2/glyc - /ox - /ldh = 0 , (1) 

(note that, to make the subsequent analysis simpler, our alge¬ 
bra differs slightly in numerical pre-factors from the one given 
in m. The overall picture is however completely equivalent). 
The glycolytic flux is in turn limited by the flux Uq that sup¬ 
plies glucose to the system, so that (7G,in < Lg, while under 
mass balance UG,in = /giyc- The crowding constraint is im¬ 
plemented as V^iyc -I- Vox + Lldh < Latp, where Vatp rep¬ 
resents the cell volume devoted to ATP production, whereas 
Lgiyc, Vox, and Vldh are the volumes taken up by the gly¬ 
colytic, OXPHOS, and LDH enzymes, respectively. Introduc¬ 
ing the constants Ogiyc, Oox, and cldh representing the vol¬ 
ume occupied per unit of ATP production by glycolytic, OX¬ 
PHOS and LDH enzymes, respectively, this can be re-cast as 

na 


tTglyc/glyc V Uqx/ox “f tiLDn/LDH ^ ^ATP , ( 2 ) 

where fl>ATP is the volume fraction available for ATP produc¬ 
tion in each cell. Using one obtains 

(olDH H- /lDH+ (uox H-/ox < ‘i’ATP ■ ( 3 ) 

We shall adhere to M in employing the empirical estimates 
Ogiyc ~ 3 X 10“^ (min/mM), oldh — 4.6 x 10“^ (min/mM), 
Qox — 2 X 10“^ (min/mM), and <i>ATP = 0.4. Equation <0 
represents the “crowding constraint”, amounting to a global 
constraint to the overall flux of metabolites the cell can invest 
in ATP production. To be consistent with na, we derived the 
crowding constraint on fluxes by limiting the volume fraction 
devoted to energy production, but the constraint 0 might as 
well be derived by limiting the amount of proteins (see, e.g., 

ESI). 

To assess the amount of ATP produced by these fluxes, we 
consider that /giyc generates 2 ATPs and 2 pyruvates per glu¬ 
cose molecule, and that fox creates 18 ATP molecules from 
each pyruvate. The conversion of pyruvate to lactate does not 
generate any ATP and it only ensures that glycolysis can con¬ 
tinue by regenerating NADH to NAD. The overall flux of ATP 
production is finally given by 

/atp = 2 /glyc + 18 /ox . ( 4 ) 

Vazquez et al’s results can be summarized as follows. Ne¬ 
glecting 0, ATP production is maximized when /giyc = Uq, 
fox = 2/glyc = 2(7g, and /ldh = 0, i.e. when all available 
glucose is taken up by the cell and used in OXPHOS. This 
solution still holds if Q is satisfied, i.e. as long as 

( 2 aox + Oglyc) Ug < <hATP ■ ( 5 ) 


If however Ug > 4* atp/ (2aox + Ogiyc) = ug, then the 
cellular resources available for high-yield energy production 
pathways has been exhausted and ATP synthesis is no longer 
limited by glucose but by intracellular resources (either vol¬ 
ume or proteins production capacity). Its maximization now 
requires that the glucose flux exceeding the possibility of pro¬ 
cessing by OXPHOS is diverted towards LDH. The maximum 
ATP production flux and the corresponding lactate flux for this 
case can be obtained by inserting /ldh = 2?7g — fox into 0 
with the constraint saturated. One gets 


and 


/ldh 


0 if Lg < MG 

aox(t^G - ^g) if Ug > UG 

“1“ ^glyc ^ ^ 
ox — ^ ^ 

^ox ^LDH 


(6) 


(7) 


2aLDH + Oglyc „ 

(TLDH = --> 0 

Uox — UlDH 


where mg = $atp/ (2aLDH + Ogiyc)- Because of mass bal¬ 
ance, /ldh also equals the lactate excretion flux so that the 
Warburg effect is described by (|^, which shows the existence 
of a threshold glucose intake above which lactate is secreted. 


Metabolic fluxes in a minimal model of two cells coupled via a 
lactate shuttle 


Let us now consider two replicas of Vazquez et al.’s cell 
and assume that the two cells can interact via the exchange 
of lactate, which, in particular, can be intaken and used as 
an alternative carbon source through reverse LDH (catalyzing 
the conversion of lactate back to pyruvate). For sakes of sim¬ 
plicity, we distinguish between a donor cell (‘don’, which can 
only excrete lactate) and an acceptor cell (‘acc’, which can 
both excrete and import lactate) (see Fig. [T]|. In order to avoid 
starvation, we furthermore impose that both cells generate a 
minimum flux /^p of ATP, i.e. 

/atp.z > /atp : * € {acc, don} (8) 

where the ATP production flux is given by 0^ namely 

/ATP.i = 2/glyc,J-f 18/ox, i /G (acc, don}. (9) 

As shown in the beginning of the previous section, when 
/giyc = Ug and G = 2(7g cells can extract the maximum 
ATP from a single glucose molecule, i.e., /atp,* = 38(7g. 
We observe that, in order to satisfy 0. the overall glucose 
supply to the system composed of two cells should satisfy 


Ug > 


Q j?min 

^Jatp 


2mg,o , 


( 10 ) 


38 
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by both glucose and lactate fluxes as /ox.acc = 2/giyc,acc + 
/LDH.don, which, vla mass balance, implies (see (|^) 

/*ATP,acc 38^glyc,acc + 18/LDH,don = 3817G,acc + 18/LDH ,don ■ 

( 12 ) 

To obtain simpler expressions, it is useful to introduce the 
shorthands 


UG,1 = UG,0 + UG , 


MG,2 


f min 

J ATP 

IScTox 


+ mg , 


(13) 

(14) 


and AUg for the amount of glucose rerouted from the accep¬ 
tor cell to the donor cell, so that 


FIG. 1. Schematic representation of the minimal model of two 
cells coupled via a lactate shuttle. Two cells (lactate donor and 
acceptor, respectively) share glucose as an energy source. Glucose 
is partitioned according to the fluxes /giyc.don and /giyc.acc. which 
convert one internal glucose molecules to two pyruvate molecules 
producing two ATP molecules. Pyruvate can then undergo oxida¬ 
tive phosphorylation (with the irreversible fluxes /ox,don and /ox,acc, 
giving 36 more ATP molecules) or LDH (with fluxes /LDH,don and 
/LDH,acc, by which lactate (L) is obtained). Lactate is for simplic¬ 
ity assumed to be secreted upon production. If both cells produce 
lactate, no coupling sets in, unless due to competition for nutrients 
under glucose limitation. If however the donor cell secretes lactate, 
the acceptor cell may intake it to replace glucose whenever its access 
to the latter is limited (e.g. because the donor cell’s glucose intake 
is large). In such cases a lactate shuttle will effectively couple the 
metabolisms of the two cells. 


which defines ug,o as the minimum glucose supply flux 
that can prevent cellular starvation. Denoting as C/G,don the 
amount of glucose available to the donor cell, if the donor cell 
maximizes ATP production, its flux organization will match 
the one found above for a single cell. It is useful to distin¬ 
guish two limiting cases: 

(a) If the donor cell employs OXPHOS exclusively, which 
happens when 2 ug,o < Ug < mg + ug,o^ have 
C^G,don = Ug - mg,o and f/G.acc = MG,o- In other 
words, the donor cell takes up all of the available glu¬ 
cose except for the amount necessary for the survival of 
the acceptor cell. In these conditions, the acceptor cell’s 
survival is guaranteed by glucose availability. 

(b) At the other extreme, if the acceptor cell can survive 
even using exclusively the lactate excreted by the donor 
cell, all glucose is intaken by the donor cell, so that 
Ug don = Ug and Ug acc = 0. This happens when 

/LDH,don > /Sp/ 18^ or 


Ug > 


f min 

/atp 

18 Qiox 


+ UG ■ 


( 11 ) 


In intermediate situations, the amount of glucose intaken by 
the acceptor cell gradually decreases from ug,o to 0 as it in¬ 
creasingly uses lactate to produce pyruvate and thus ATP. In 
such cases, the oxidative pathway of the acceptor cell is fed 


C^G,acc = MG,0 - AUg 
C^G,don = Ug — MG,0 + AUg ■ 

Substituting UG,a.cc and imposing minimal ATP requirement 
in ( [T 2 I 1 , we obtain 

00 

/LDH,don = 

Since the donor cell maximizes ATP production, we obtain 
a second expression for /LDH.don as a function of AUg by 
substituting [/G,don: in 

/LDH,don = aox(C^G,don —Mg) = CtoxiUG—UGfi + AUG — UG) ■ 

This, combined with ( [Tbl l, determines the extra amount of glu¬ 
cose available to the donor cell 

38/18 

Expression is valid if mg,i < Ug < mg, 2 - For Ug < 
UG,i one has AUg = 0, while when Ug > mg ,2 AUg = 
mg,o- We can therefore summarize the fluxes for the donor 
cell as 


fo 

if 2mg,o <Ug< mg,i 

/LDH,don = < aox(C^G + AUg “ 

mg,i) if mg,i < 17g < mg,2 

[aox(C^G — Mg) 

if Ug > Mg.2 


(18) 

and 


f 2(17g - mg.o) 

if 2ug.o <Ug< mg,i 

fox,don = < aLDH(MG — t^G,don) 

if mg,i <Ug < Mg.2 

[aLDH(MG - Ug) 

if Ug > mg,2 


(19) 

Because the donor cell uses up most of the glucose, the ac¬ 
ceptor cell stays at the minimum ATP production rate /^p 
for Ug < Mg, 2 - In this regime, the fluxes of the acceptor cell 
are fixed and can be derived from the ones of the donor cell. 
One gets 


/LDH,acc 


0 if 2 ug,0 <Ug< Mg,1 

—/LDH,don if Mg,1 < Ug < Mg,2 


( 20 ) 
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FIG. 2. Qualitative behaviour of the single-cell and of the donor- 
acceptor system for the coarse-grained model, (a) In the cell- 
autonomous model, lactate overflow occurs when the glucose intake 
overcomes a threshold. Correspondingly, the flux through oxidative 
metabolism increases until the threshold before slowly decreasing 
once the crowding constraint is saturated and fermentation sets in. 
(b) In the donor/acceptor system, the donor behaves essentially as an 
autonomous cell and the acceptor adapts to it. For low glucose in¬ 
takes, it operates oxidative pathways at small rate. When the donor 
saturates the crowding constraint excreting lactate, the acceptor im¬ 
ports it and uses it at a substrate to increase the flux through oxidative 
metabolism. 


and 


Metabolic fluxes and crowding constraint in a cell-autonomous 
large scale metabolic network 

In order to characterize the exchange of carbon equivalents 
among interacting cells more in detail, we have analyzed a 
human metabolic network derived from the reconstructed hu¬ 
man reactome ED, to which we refer as the ‘Human Core 
Catabolic Network’ (HCCN, see ‘Materials and Methods’). 
The HCCN is composed by 67 metabolites and 75 reactions 
(including uptake fluxes), 23 of which are reversible. The 
crowding constraint that accounts for finite cellular resources 
is represented as 

«glyc/HEX + aox(/pDH + /gLUn) + OLDh/lDH < ‘J’ATP 

(25) 

where /hex denotes the Hexokinase-1 flux, which irre¬ 
versibly channels glucose into glycolysis, /pdh denotes the 
flux through pyruvate dehydrogenase, by which pyruvate is 
diverted into the TCA cycle and /glun denotes the mithocon- 
drial flux of nitrogen metabolism through glutaminase. We 
explore the space of flux patterns consistent with the con¬ 
straints by sampling solutions of the mass balance equations 
Sf = 0 (S denoting the stoichiometric matrix and f being a 
flux vector) such that flux vectors fdon of the donor cell ap¬ 
pear with probability 


f _ / 2'“g,o if 2wg,o <Ug<ug 

[2{7G,acc +/LDH,don if Mq.I < Ug < UG,: 


For large total glucose intakes, the donor cell excretes more 
lactate than necessary for bare survival of the acceptor cell, 
which can produce any amount of ATP as long as it is larger 
than /^p. The metabolic state of the acceptor cell is there¬ 
fore not uniquely defined. When Ug > ug, 2 , for the fluxes of 
the acceptor cell one finds 


—/pDH.don < yLDH,acc ^ 

/ox,acc — ^LDH,acc ; 

while the ATP production is given by 

/atp ,acc = —18/LDH,acc 


•min 

ATP 


38 


( 22 ) 

(23) 

(24) 


Note that /pDH.acc < 0 because the acceptor cell intakes lac¬ 
tate and reverses LDH. 

The solution to this model is described pictorially in Fig. 
and discussed in Fig.[^ where we also show the feasible range 
of values of the fluxes of the acceptor cell. 

In essence, this coarse-grained model suggests that an im¬ 
balance in the energetic demands of two cells can induce a 
metabolic coupling driven by a lactate shuttle from the high- 
to the low-demand cell, and provides a qualitative scenario 
for how carbon utilization in both cells changes by tuning the 
overall glucose supply. We shall now see that such a pic¬ 
ture is fully recovered within a large-scale model of human 
metabolism. 


-P(fdon) = A explPfATP ,don] ; (26) 

where /3 > 0 is a parameter and A is a normalization con¬ 
stant. The rationale is that, by tuning the value of /3, we can 
pass from an unbiased sampling in which /atp .don takes on 
every allowed value with uniform probabilty (P = 0) to one 
in which the donor cell maximizes its ATP output (/3 —oo), 
thereby obtaining a complete and refined picture of how dif¬ 
ferent degrees of optimization by the donor cell impact the 
flux pattern of the acceptor cell. Solutions at the two extremes 
can be obtained with standard methods, as linear program¬ 
ming, but only our method allows us to sample realistic sub- 
optimal optimizations. In Fig. S6 we show how the ATP pro¬ 
duction increases as a function of /3. In the following, as a 
prototype for sub-optimal optimization, we choose /3 = 5 that 
corresponds to roughly 70% of the maximum ATP production. 
Model and algorithmic details are discussed in ‘Materials and 
Methods’. 

In Fig. 1^ we show the ATP production of a single HCCN 
cell as a function of the available glucose. If ATP produc¬ 
tion is maximized (large /3), one first encounters a regime 
with a yield of roughly 30 ATPs per glucose molecule (to be 
compared with the theoretical value of 38) and then a regime 
where the yield is about 1.7 ATPs per glucose. ATP is pro¬ 
duced by OXPHOS in the former case, and by fermentation 
in the latter (see also Fig. and Fig. |^). As expected, net¬ 
works maximizing ATP shift their metabolic strategy at the 
glucose intake for which the resources available for ATP pro¬ 
duction by oxidation is exhausted. Note that for /3 = 5 one 
obtains an ATP efficiency very close to optimal. Quite im¬ 
portantly and surprisingly, however, for /3 = 0 (i.e. when 
no ATP maximization is performed) the emerging picture is 
qualitatively preserved, albeit with lower ATP yields. Indeed 
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Ug 


(b) ATP production 




Ug 


^ATP.don 


(c) Fermentation and OXPHOS fluxes 


(d) Metabolic switch in donor cell 


FIG. 3. Solution of the minimal model of two cells coupled via a lactate shuttle. A lactate donor cell that maximizes ATP production is 
coupled to an acceptor cell, (a) Glucose intake for donor and acceptor as a function of total glucose supply, (b) ATP production as a function of 
the total glucose supply. The shaded area indicates that all values within that region are feasible, (c) Flux through fermentation /ldh (circles) 
and oxidative phosphorylation /ox (crosses) as a function of the total glucose supply, (d) Fraction of ATP produced via fermentation (black) 
or via oxidative phosphorylation (red) in the donor cell as a function of the total ATP it produces. 


solutions sampled at /3 = 0 appear to employ a mixture of 
OXPHOS and fermentation even at low glucose intakes. It is 
remarkable that the resources-driven shift still occurs at uq, 
as in this case it corresponds to a largely suboptimal value 
of ATP production. In other words, the HCCN might devote 
cellular resources to increase ATP production, but to do so 
it must explicitly be pursuing ATP maximization. This im¬ 
plies that the crossover from a high- to a low-yield strategy is 
a robust, embedded property of the network (and of the con¬ 
straints imposed) and not an exclusive characteristic of the ex¬ 
tremal solution that maximizes the ATP production. 

In Fig.|^, we focus on the fluxes through PDHm and LDH, 
that indicate whether pyruvate is directed towards OXPHOS 
or fermentation. We see that high ATP yields are obtained 
by using OXPHOS exclusively. The typical network state at 
/? = 0 is less efficient at all glucose intakes, as it always di¬ 
verts a larger fraction of pyruvate to fermentation compared 
to the ATP-maximizing cell. Finally, in Fig.|^, we detail how 
glucose partitions between fermentation and OXPHOS. The 
network displays a reversal of glucose fate as the crowding 
constraint is saturated independently of whether ATP produc¬ 
tion is being maximized or not. The latter case however turns 
out to be genetically less efficient. 


Metabolic fluxes in a large scale metabolic network of two cells 
coupled via a lactate shuttle 


We now consider two replicas of the HCCN and again dis¬ 
tinguish between a lactate acceptor and a donor cell. We also 
impose that there is no external lactate source, and explore 
the scenario where the donor optimizes ATP production to a 
degree tuned by the parameter (3. Fig. shows that, as ex¬ 
pected, an ATP-maximizing donor cell sequesters all of the 
available glucose except for the small amount required for the 
acceptor’s survival. This is reflected by the ATP production 
curves as a function of the total glucose supply displayed in 
Fig. 1^. In such conditions, the acceptor’s ATP production 
fluxes matches the minimum required for survival, i.e. /^p 
(which is set to be equal to 1 uq for simplicity) until the donor 
switches to aerobic glycolysis, thereby excreting lactate. The 
donor, on the other hand, goes through an initial phase of ex¬ 
clusive OXPHOS use, followed by a switch to aerobic glycol¬ 
ysis when the crowding constraint does not allow for a further 
increment of the mitochondrial flux. In Fig. one indeed 
sees that glucose is mainly channeled to OXPHOS as long as 
the crowding constraint allows for it. As soon as the latter is 
saturated, the donor diverts pyruvate to LDH, usefully regen¬ 
erating NAD from NADH and thus avoiding glycolysis halt. 
LDH generates lactate, which is then expelled and intaken by 
the acceptor. As shown in Fig. |^, this is accompanied by a 
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(b) Fermentation and OXPHOS fluxes 



FIG. 4. An isolated HCCN maximizes the ATP yield by directing 
glucose to OXPHOS even in absence of an active ATP maximiza¬ 
tion. (a) Average ATP production as a function of the re-scaled av¬ 
erage glucose supply, (b) Average flux through LDH (crosses) and 
PDHm (circles) as a function of the re-scaled average glucose sup¬ 
ply. (c) Fraction of ATP produced via glycolysis (crosses) or via 
OXPHOS (circles) as a function of the re-scaled total ATP produced. 
The flux through each pathway is re-scaled by half the amount of glu¬ 
cose intaken by the cell (because with one molecule of glucose cells 
produce two molecules of pyruvate). Curves describe the behaviour 
of an ATP-maximizing HCCN (black lines, /3 = 50), a loosely max¬ 
imizing donor (red lines, /3 = 5) or the result of a uniform sampling 
of the allowed flux states for a HCCN (blue lines, /3 = 0). Error bars 
for the s.e.m. are smaller than symbols. 




FIG. 5. When a lactate donor maximizes its ATP production it 
intakes most of the glucose supplied to the two-cell system. The 
ATP production by the acceptor cell increases in correspondence 
to a decrease in efficiency of the donor’s metabolism, (a) Glu¬ 
cose intakes for two coupled cells as a function of the total glucose 
available to the donor-acceptor pair, (b) ATP produced by the donor. 
Curves describe the behaviour obtained for two coupled HCCN cells 
with an ATP-maximizing donor (black lines, j3 = 50), a loosely 
maximizing donor (red lines, /3 = 5) or for an unbiased sampling of 
the two-cell solution space (blue lines, /3 = 0). Error bars for the 
s.e.m. are smaller than symbols. 


reversal of LDH in the acceptor; lactate is transformed into 
pyruvate that is then channeled to OXPHOS through PDHm. 
The acceptor can thus spare some cellular resources to pro¬ 
duce pyruvate, at the cost of becoming strongly dependent on 
the donor for ATP production. 

If one instead averages over all solutions without biasing 
for ATP production by setting /3 = 0, the two cells produce 
comparable amounts of ATP and share glucose more evenly 
(see Fig. |^). The fact that the donor turns out to employ 
slightly more glucose than the acceptor is due to the intrin¬ 
sic asymmetry that is introduced by not letting the donor in¬ 
take lactate. This asymmetry is also sufficient to drive, even 
when an unbiased statistical picture of the solution space is 
considered, a net lactate exchange as the typical behaviour 
Both cells, however, produce a factor 6 (roughly) less ATP 
than the maximum possible, given the glucose influx. Despite 
a similar overall ATP production profile, the internal fluxes 
of donor and acceptor differ significantly. While the donor 
distributes glucose almost evenly through LDH and PDHm, 
the acceptor can only generate ATP via OXPHOS, since it 
intakes lactate as an extra source of carbon. In Supplemen¬ 
tary Fig. SI, we show that two symmetric cells with unbiased 
ATP production show identical glucose intake, ATP produc¬ 
tion, and internal flux distribution. Even for symmetric cells, 
however, when the donor optimizes the ATP production, we 
observe lactate shuttling between donor and acceptor cells, as 
illustrated again by Supplementary Fig. SI. In summary, for 
maximal ATP production (/3 —oo) there is no difference be¬ 
tween ATP production and internal flux distributions of sym¬ 
metric and asymmetric cells-couples. For simplicity in the 
analysis of the data, we used the asymmetric case where the 
donor cell cannot intake lactate. 

We observe that the fraction of ATP produced via oxidative 
phosphorylation utilized by the acceptor cell differs substan¬ 
tially from the one of the donor cell. In Fig. we display 
the relative usage of the fermentative and the oxidative path¬ 
way and observe that, for low glucose supply, the acceptor cell 
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(a) Fluxes in acceptor cell 




(c) Carbon partitioning in donor cell (d) Carbon partitioning in acceptor cell 


FIG. 6. Oxidative and fermentative fluxes in a donor-acceptor system, (a)-(b) Average flux through LDH (circles) and PDHm (crosses) as 
a function of the average glucose supplied to the two-cell system, (c)-(d) Fraction of ATP produced via glycolysis (circles) or via oxidative 
phosphorylation (crosses) as a function of the total ATP produced hy the donor cell. Pathways are normalized by the “pyruvate-equivalent”, i.e. 
by the sum of lactate and half-glucose intakes. Curves describe the behaviour obtained for two coupled HCCN cells with an ATP-maximizing 
donor (black lines, P = 50), a loosely maximizing donor (red lines, /? = 5) or for an unbiased sampling of the two-cell solution space (blue 
lines, p = 0). Error bars, which represents s.e.m., are smaller than symbol sizes. 


only employs the oxidative pathway, independent of whether 
the donor cell maximizes ATP or not. When the donor cell 
maximizes ATP production, it sequesters most glucose and 
the acceptor cell is therefore forced to feed on lactate, which 
can only be usefully converted to ATP by means of OXPHOS. 
Conversely, when the donor cell does not maximize ATP pro¬ 
duction, it secretes a sizable amount of lactate because its 
metabolism is inefficient. The lactate acceptor, however, is 
also mainly feeding on glucose (see Fig.|^) and could in prin¬ 
ciple secrete lactate just like the donor cell. The reason why a 
purely oxidative metabolism is observed lies in the reversibil¬ 
ity of LDH. Although lactate intake by the acceptor is small, 
it suffices to force the pyruvate flux towards TCA and OX¬ 
PHOS, thereby making the acceptor more efficient than the 
donor in energetic terms. 

In order to assess the robustness of the occurrence of lactate 
overflow metabolism under the imposed constraints, we have 
further analyzed the flux configurations in a dono/acceptor 
system in which the donor maximizes the production of 
biomass precursors. We show here (see Fig. 0 results ob¬ 
tained by maximizing the biomass defined as in Table S9, 
which essentially reproduce those described above. Such a 
robustness becomes less surprising in the light of the fact that 
the qualitative features of the switch from oxidative to non- 
oxidative metabolism are obtained even in an unbiased sam- 



FIG. 7. Donor’s biomass flux and lactate overflow, and accep¬ 
tor’s lactate intake for two coupled HCCNs with a biomass- 
maximizing donor. The qualitative behavior obtained when the 
donor is maximizing the ATP flux is reproduced in a more realis¬ 
tic case in which a biomass objective function is considered. See 
Table S9 for the biomass coefficients. 

pling of the steady states. This observation suggets that such 
a scenario is to a large degree embedded in the stoichiometry 
and in the main topological features of the underlying reaction 
network. 
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Correlation coefficients among metabolic fluxes of two cells 
coupled by lactate shuttle under glucose limitation 

To formally assess the extent of metabolic coupling be¬ 
tween an ATP-maximizing lactate donor and a lactate accep¬ 
tor, we computed the matrix of normalized Pearson correla¬ 
tion coefficients of each pair of fluxes in the solutions sampled 
for different levels Uq of the glucose supplied to the system. 
The Pearson coefficient between random variables X and Y 
is defined as r = cov(X,Y)/((7xC’'y), where cov{X,Y) de¬ 
notes their covariance and ax and ay stand for their respec¬ 
tive standard deviations, r ranges from —1 to 1 and quantifies 
the linear dependence of the two variables. More precisely, 
the linear correlation between variables X and Y is more pos¬ 
itive the closer r is to 1, and more negative the closer r is to 
— 1, while X and Y can be considered uncorrelated if r ~ 0. 
For sakes of clarity, we have considered the correlations aris¬ 
ing in a system formed by an ATP-maximizing lactate donor 
(large /3) and a lactate acceptor. For smaller values of /3 corre¬ 
lations get weaker while maintaing the same qualitative struc¬ 
ture. Fig. [^displays three reduced correlation matrices (ob¬ 
tained for three different values of the overall glucose supply) 
where only a small subset of fluxes (each presentative of a 
different biologically relevant pathway) appears. Full matri¬ 
ces for three choices of Uq are instead shown in Figs. S3-S5. 

Intuitively, correlation patterns should depend strongly on 
Uq. The essence of this dependence, which clearly emerges 
from Fig. is that cross-correlations between donor and ac¬ 
ceptor build-up as the glucose supply increases, are maximal 
when the acceptor’s energetics depends on donor-derived lac¬ 
tate, and then decrease again rapidly to zero when the glucose 
supply is large enough to sustain both cells. 

If Uq is low (sub-threshold for lactate overflow), the 
ATP-maximizing donor necessarily outcompetes the accep¬ 
tor for the available glucose, establishing a degree of cross¬ 
correlations that is driven mainly by the unbalanced nutrient 
partitioning. The donor indeed employs OXPHOS almost ex¬ 
clusively, while still leaking out a small amount of lactate 
which is taken up by the acceptor (i.e. there is no accumu¬ 
lation in the tissue). This suffices to establish the overall cor¬ 
relation among cells that can be seen in Fig. [^. 

As the glucose supply is increased further (above thresh¬ 
old for lactate overflow), the ATP-maximizing donor switches 
to aerobic glycolysis and secretes lactate that is shuttled to 
the acceptor in large amounts (viz. the strong negative cor¬ 
relation arising between the donor’s lactate outflux and the 
acceptor’s lactate influx). In such conditions, the overall 
cross-correlations increase, see Fig. [^. In particular, gly¬ 
colytic fluxes in donor and acceptor anti-correlate, oxidative 
pathways weakly anti-correlate, while te donor’s biosynthetic 
pathways correlate weakly with the acceptor’s oxidative path¬ 
ways. Notice also that ATP flux in the donor weakly anticor¬ 
relates with every pathway in the acceptor’s metabolism. 

At high enough Uq, finally, cross-correlations decay as 
both cells can rely entirely on glucose for their energetics. 
Because the maximum glucose uptake for an HCCN is given 
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(a) Low glucose 
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(b) Intermediate glucose 
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(c) High glucose 


FIG. 8. Restricted matrices of Pearson correlation coefficients 
for two coupled HCCN cells when the lactate donor maximizes 
the ATP production and the overall glucose supply is large shows 
that the two cells are not correlated. The intensity of the color rep¬ 
resents the magnitude of the correlation coefficient (see scale on the 
right hand side). The two cells can independently access glucose and 
internal fluxes of the lactate acceptor and donor are essentially un¬ 
correlated. The four representative reactions displayed for each cell 
are the glucose influx (Glc, a proxy for glycolytic activity), LDH (a 
proxy for lactate overflow and exchange), PDH (a proxy for oxidative 
metabolism) and the ATP production flux (ATP). 






10 


by = $ATp/ (qg + Scldh) (this is determined by the 

crowding constraint ( |25| l when glucose is entirely channeled 
to fermentation, i.e. when /hex = (^G,don = /ldh/2), one 
may assume that Ug is “large enough” when it is larger than 
2(7 q as both cells would in this case have access to as much 
glucose as they can intake. In this case, no donor-acceptor 
correlations should be expected. Fig. shows that this is in¬ 
deed the case: intra-cellular correlations, described by the di¬ 
agonal blocks, greatly exceed inter-cellular ones (off-diagonal 
blocks). 

In summary, when glucose availability is limited, the fast¬ 
growing donor necessarily outcompetes the acceptor for glu¬ 
cose and a network of cross-correlations is established be¬ 
tween the metabolic systems of the two cells. Metabolic path¬ 
ways in donor and acceptor correlated however most strongly 
when a lactate exchange from donor to acceptor is established, 
and decay again as the two cells become effectively decoupled 
in glucose-rich media. The build-up of a metabolic partner¬ 
ship is perhaps best seen by the fact that the donor’s glucose 
influx C/G,don develops a correlation with the acceptor’s LDH 
flux at intermediate values of Uq. Note also that C/G,don cor¬ 
relates more and more strongly with the donor’s LDH flux as 
Uq increases. 

DISCUSSION 

Experimental evidence on lactate shuttles in cancer 

'While lactate shuttling is a likely scenario in many physi¬ 
ological conditions, from muscle cells undergoing intense ac¬ 
tivity II22II to neuron-glia energetics Il23l426ll . and has been re¬ 
ported to take place intra-cellularly at the mitochondrial and 
peroxisomal membrane i22\ . no direct measurement of inter¬ 
cellular lactate exchange in cancer exists. Strong indirect evi¬ 
dence, however, suggests that such a scenario is indeed plau¬ 
sible. 

Indeed, tumor-stroma and tumor-tumor metabolic inter¬ 
actions are currently being characterized at various levels. 
The emerging picture increasingly suggests that cancer sus¬ 
tainment is a non-cell-autonomous phenomenon Il27l430l and 
that stromal cells might be potential targets for cancer treat¬ 
ment ED. Lactate exchange in particular has been investi¬ 
gated both as fueling oxidative cancer cells and as support¬ 
ing non-aberrant cancer-associated cells ll32l[3^ . The latter 
case corresponds to the model discussed here. Signatures of 
tumor-to-stroma lactate shuttle have been reported in terms of 
the over-expression of monocarboxylate lactate transporters 
jointly with increased PDH activity in cancer-associated fi¬ 
broblasts (CAFs) lfT9l[34l. Likewise, tumor-derived lactate 
has been found to play a major signaling role (specifically, 
for the initiation of tumor angiogenesis) in vascular endothe¬ 
lial cells II35I . Inhibition of the lactate transporter MCTl has 
therefore been proposed as a possible anti-angiogenic strategy 

m- 

Taken together, these observations suggest that tumor cells 
and their associated stromal and vessel cells can be seen as 
a collective, synergistic metabolic unit where each compart¬ 
ment carries out complementary functions reflected in their 


energetic strategies. In such a microenvironment, aerobic stro¬ 
mal cells serve an essential ‘bioremediative’ role by remov¬ 
ing potentially toxic compounds, thereby reducing acidity and 
positively feeding back on anaerobic cancer growth. The ex¬ 
istence of a consistent imbalance in energetic demands across 
different compartments is crucial to establish this scenario. 

Intercellular shuttling of lactate towards oxidative cancer 
cells is known to occur in two distinct forms, namely either 
from non-oxidative (e.g. hypoxic) to oxidative tumor cells 
ll3^ or through the ‘reverse Warburg effect’, i.e. the onset 
of aerobic glycolysis in CAFs IIJTHtTII . The former case de¬ 
scribes the metabolic sysmbiosis that is established e.g. be¬ 
tween more and less hypoxic regions of a tumor, which puta¬ 
tively allows glucose to be delivered more efficiently to more 
hypoxic regions 14^ (see ll43l for a game-theoretic approach 
to modeling this kind of tumor-tumor interplay). In the lat¬ 
ter scenario, cancer cells induce oxidative stress in CAFs, re¬ 
sulting in CAFs switching to anaerobic glycolysis Il44l . The 
secreted lactate is then imported by cancer cells for use in aer¬ 
obic pathways. This effect has indeed been observed in vitro 
EH HD, and its relevance in vivo is currently under scrutiny 
ll3^ . Clearly, however, this type of shuttling is not necessarily 
triggered by a strong imbalance in energetic demands between 
the lactate donor and acceptor cells, and different constraints 
in metabolic activity are likely to be required to understand its 
origin within genome-scale models. 

It is worth noting that recently developed microfluidic plat¬ 
forms allow to probe the tumor microenvironment with un¬ 
precedented resolution and control over nutrient supply ll45l . 
More light will hopefully be shed on its activity as a functional 
metabolic assembly and on the role of lactate in specific. 

Going beyond cancer, however, such studies highlight the 
fact that replacing the natural environment in which cells live 
with a cell culture that is possibly optimized for the cell’s 
needs might severely limit our ability to understand the ac¬ 
tual behavior of the system in vivo. On the other hand, they 
suggest that cell-autonomous models may be unable to cap¬ 
ture some essential features of the energetics of cells: when 
cells employing different energetic demands share a limited 
resource, a mutually beneficial microenvironment character¬ 
ized by a large-scale exchange of chemical species can be es¬ 
tablished. A thorough understanding of cellular growth strate¬ 
gies should take these aspects into account. 


Conclusions 

In this article we have studied the lactate-driven coupling 
that is established between cells with different energetic re¬ 
quirements, showing that a lactate shuttle appears robustly 
as a consequence of basic physico-chemical constraints. Our 
model covers time scales over which the cellular metabolic 
networks have to adapt to the establishment of an imbalance 
in energy demands (as well as, for instance, in the levels of 
main carbon transporters), which are much shorter than the 
time required to alter the nutrient availability profile signifi¬ 
cantly (e.g. via vascularization). We have illustrated the emer¬ 
gent metabolic interaction scenario in a simplified model that 
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captures the essential fact that, for two such cells, recycling 
fermentation products may be a good strategy to optimize cel¬ 
lular glucose usage. To reach these conclusions, we have con¬ 
sidered a human metabolic network model of cells interacting 
via lactate shuttle and sampled its feasible metabolic states ac¬ 
cording to a prescribed probability distribution by which we 
can tune the degree to which one of the two cells is pushing 
its ATP production. The scenario obtained in a highly sim¬ 
plified (but exactly solvable) model is fully retrieved in this 
case, where a thorough analysis of correlation patterns reveals 
further details of the cell-to-cell coupling. 

Such a coupling would have serious physiological implica¬ 
tions. First, lactate recycling implies that lactate may accu¬ 
mulate in tissues only in late-stage tumors. Before neighbor¬ 
ing cells saturate their processing capabilities, the lactate ex¬ 
pelled by cancer cells would be taken in by non-aberrant cells, 
suggesting that localized lactate measurements might identify 
early-stage tumors. Second, lactate shuttling and recycling 
systems should be regarded as potential targets for treatment. 

There are many aspects of the approach presented here 
that can be further dissected. First is the nature of con¬ 
straints. From a physical point of view, imposing “volume” 
constraints of the type considered here is equivalent to impos¬ 
ing finite flux capacities on energy-efficient pathways. There 
are however strong indications that unicellular organisms ac¬ 
tively down-regulate energy-efficient pathways at high nutri¬ 
ent levels iniia. The fact that cells invest energy and re¬ 
sources in silencing the synthesis of enzymes involved in aer¬ 
obic pathways suggests that the switch to an energy-inefficient 
pathway provides a real, physiological advantage for the cell. 
Its origin and nature are currently not understood, although the 
inclusion of additional “costs” due to regulation in metabolic 
models strongly suggests that the observed phenomenology 
can be captured, at least to some degree, by accounting for the 
fact that a shift in energetic strategy requires a sizeable change 
of a cell’s protein repertoire m. 

It is also important to note that the ETC in eukaryotes 
takes place in a separate cellular compartment (the mitochon¬ 
drion), at odds with prokaryotes like bacteria, where it oc¬ 
curs is the cell’s periplasm. We have entirely neglected the 
complications due to spatial organization, which is reasonable 
as long as one wants to focus on the ATP yield of energy- 
producing strategies. However, this aspect will become im¬ 
portant for higher-resolution modeling at genome scale. Like¬ 
wise, a more refined spatial modeling will allow for an in sil- 
ico analysis of a possible role of intra-cellular lactate shut¬ 
tling in cancer ll22]l . In brief, the key idea behind this is that 
glycolysis-derived lactate could be employed directly as an 
additional energy source for mitochondrial oxidation instead 
of being excreted. Elementary biochemical properties of the 
LDH reaction indeed would appear to favour lactate (rather 
than pyruvate) production in the cytosol. Moreover, experi¬ 
mental facts indicate a possible role of intra-cellular lactate 
shuttling in cardiomyocites, neurons and skeletal muscle cells 
facing high energy demands HSi . Accounting for it in a com¬ 
putational scheme may help reveal novel details about possi¬ 
ble pathways of energy production during physical exercise or 
in rapidly growing cells. 


Eurtermore, while this work has focused on lactate as a 
key mediator of metabolic, energy-driven intercellular inter¬ 
actions, it is worth stressing that recent studies focusing on 
the comprehensive analysis of the exometabolome of growing 
cells has revealed that metabolic interactions encompass many 
more chemical species and is tightly coupled with the growth 
regime. In bacteria, for instance, released intermediates in¬ 
clude not only carbon equivalents (like acetate, pyruvate and 
ethanol) but also amino acids and central metabolic interme¬ 
diates [like fructose-6-phosphate, glucose-6-phosphate, 2/3- 
phosphoglycerates, phosphoenolpyruvate, acetyl-CoA, citrate 
and a-ketoglutarate], higher amounts of which have been 
found to be consistently released when carbon availability is 
high, while intermediates are typically in-taken in carbon- 
limitation El. Similar interactions are now known to occur 
in cancer: leukaemia cells have indeed been recently shown 
to employ cysteine derived from bone-marrow stromal cells 
as a means to fight oxidative stress ll^ . As the scenario un¬ 
derpinning the establishment of such couplings is increasingly 
elucidated, it will become possible to set up more refined and 
realistic models to capture a greater extent of their physiolog¬ 
ical relevance. 

Einally, it should be kept in mind that aerobic glycolysis 
in cancer cells, especially at later stages of oncogenic devel¬ 
opment, could be due to the fact that the tumor microenvi¬ 
ronment is hypoxic. The model we consider here addresses 
a faster time scale, over which the energy balance of neigh¬ 
bouring cells disrupts while the overall amount of resources 
available remains (roughly) unchanged, and suggests that at 
these stages tumor-to-stroma metabolic interactions can pro¬ 
vide a mutually beneficial solution to the imbalance. It is also 
worth mentioning that cases are known in which malignant 
cells rely on increased oxidative phosphorylation rather than 
on aerobic glycolysis for energy production (as found, for in¬ 
stance, in transformed human mesenchymal stem cells 151). 
Taken together, these results suggest that cancer’s bioenerget¬ 
ics, starting with the aerobic glycolysis versus oxidative phos¬ 
phorylation “dilemma”, may be a dynamically modulated pro¬ 
cess that differs widely across tumor types ^491 

The modeling approach discussed here, based on explor¬ 
ing the solution space induced by mass-balance equations 
rather than on optimizing a prescribed objective function, is 
in our view the most suited to deal with multi-cell systems 
in which extended cell-to-cell metabolic interactions are es¬ 
tablished by one cell’s deregulated metabolic demands. The 
sampling method employed here is scalable and easy to im¬ 
plement, providing a highly promising tool for further stud¬ 
ies. It would in particular be important to analyze the emer¬ 
gence of these effects in full-fledged genome-scale models of 
specific cancers, as can be obtained e.g. by refining recon¬ 
structions of human metabolic networks with cancer expres¬ 
sion data. The full complexity of metabolic cell-to-cell in¬ 
teractions has only just started to be uncovered. Understand¬ 
ing its functional relevance by detailed in silica models may 
provide new insight into the mechanism of cancer progression 
and further advance the search of specific, targeted treatments. 
It would be important to assess the relevance of lactate shut¬ 
tling in the context of tumors through experiments probing 
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metabolism (or metabolite exchanges) directly, rather than by 
collecting indirect evidence in the up- or down-regulation of 
specific transporters, and metabolic solution space sampling 
is known to provide useful keys for experiment design ll53l . 
Our results indicate at least two potential ways to obtain use¬ 
ful information. One possibility requires setting up an assay 
where the ability of non-aberrant cells to intake lactate can 
be modulated in the presence of bona-fide lactate-secreting 
cancer cells. Such a modulation could be achieved, for in¬ 
stance, by either varying the relative concentration of donor- 
and acceptor-cells, or by changing the expression of lactate 
and glucose transporters in supporting cells. It is reasonable to 
expect that, in such a setup, the extracellular acidification rate 
due to the accumulation of lactate in the extracellular medium 
should negatively correlate with the lactate removal capacity 
of the supporting cells. On the other hand, results from our 
statistical sampling of the solution space suggest that differ¬ 
ent pathways (and hence, likely, the expression levels of dif¬ 
ferent enzymes) should be tuned in a specific way by chang¬ 
ing the overall glucose supply in a mixed culture shared by 
aberrant and non-aberrant cells. Monitoring expression lev¬ 
els in different, controlled nutrient conditions would therefore 
provide key validation to the coupling picture discussed here 
(and, more generally, to any cell-to-cell coupling scenario for 
cancer sustainment). 


MATERIALS AND METHODS 
Human core catabolic network for a single cell 

We built a realistic metabolic network of ATP production 
from the decompartmentalized human reactome Recon-1 ED 
that we use throughout the manuscript. (While the more re¬ 
cent Recon-2 reactome ||50l provides a higher resolution re¬ 
construction of human metabolism, the degree of detail pro¬ 
vided by Recon-1 more than suffices for our purposes.) The 
HCCN includes the following Recon-1 pathways: Glycolysis, 
Pentose Phosphate Pathway, Citric Acid Cycle, and Oxidative 
Phosphorylation. In addition, we included anabolic pathways 
for the production of biomass precursors (amino acids and 
fatty acids) in an effective way. Details of the reconstructed 
network are given in Supporting Text SI. We then proceeded 
by removing the leaves of the resulting network (a necessary 
pre-processing step required to implement the mass-balance 
constraints described below). To rid the model of a priori in¬ 
feasible loops, we resorted to the method described in ED. 
A single infeasible cycle was identified, which was fixed by 
leaving only one of the two isoforms of the isocitrate dehy¬ 
drogenase that uses NADP as cofactor. 

The final HCCN we employed in this study is composed 
of 67 chemical species (listed in Table S2) and 75 reactions 
(listed in Tables S3 and S4), including the uptake/secretion 
of 11 metabolites (namely molecular oxygen, carbon dioxide, 
water, hydrogen, lactate, glucose, ammonia, phenylalanine, 
methionine, glutamine, and methyl group—which represents 
a generic methylation-whose bounds are fixed as described 
in Table S5), and reactions consuming ATP, aminoacids and 


palmitic acid. Among intracellular reactions, 23 are reversible 
according to Recon-1 ’s thermodynamic assignments. We have 
considered a medium with variable maximum glucose in¬ 
take, fixed maximum glutamine intake and unbounded oxygen 
availability. Finally, ATP consumption presents a lower bound 
standing for the minimum ATP consumption flux /^p nec¬ 
essary for cell survival (the bounds of all reactions are listed 
in Table S6 and S7). 


Human core catabolic network for two coupled cells 

In order to analyze the coupling of two identical cells, we 
simply replicated the HCCN twice, the only difference ly¬ 
ing in the uptake reaction for lactate. In specific, the lac¬ 
tate donor can only secrete lactate (through an irreversible 
reaction), while the acceptor can both excrete and import 
it. The coupling is ultimately determined by the partition¬ 
ing of glucose and by the shuttling of lactate. The former 
is a shared resource, while the latter can be exchanged be¬ 
tween the two cells: both cells can produce lactate, but there 
is no external source of lactate. We have included two fur¬ 
ther constraints for the sum of glucose and lactate flux of 
the donor-acceptor system, i.e., f/c.don + t^c.acc < Uq and 
C^LAC.don + t/LAC.acc > 0 (see Table S8). As for a single 
HCCN, glucose can only be imported while lactate can only 
be excreted. 

Crowding constraint for the human core catabolic network 

The crowding constraint within the HCCN has been imple¬ 
mented as 

Oglyc/HEX + aox(/pDH + /GLUN) + OLDh/lDH < ^’ATP 

(27) 

where /hex denotes the Hexokinase-1 flux, /pdh denotes 
the flux through pyruvate dehydrogenase, by which pyruvate 
is diverted into the TCA cycle, and /glun denotes the flux 
through glutaminase. For Ogiyc, oldh and apDH we use the 
numerical values given in m for a coarse-grained model, 
since for unitary flux both /hex and /ox transport the same 
amount of carbons present in glucose, which is exactly what 
the effective fluxes defined in llT4l do. 

Inequality ( |27] i is valid for a cell that excretes lactate. When 
lactate can also be intaken, as occurs for the coupled HCCN 
cells, the flux through LDH can become negative. The more 
general form of the crowding constraint we consider is there¬ 
fore 

Oglyc/nEX + aox(/pDH + /glun) + OLDhI/lDhI < ‘hATP. 

(28) 

For both the single HCCN and the two-HCCN models, we 
provide as supplementary materials files containing the lists 
of metabolites and reactions (with bounds), and the explicit 
expressions of each reaction in the HCCN in terms of inde¬ 
pendent variables, which provide a full characterization of 
the solution space polytope. The files can also be down¬ 
loaded from http://chimera.romal.infn.it/SYSBIO where an 
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SBML file with the model and a text file with the stoichio¬ 
metric matrix of the single cell HCCN are also made avail¬ 
able. HCCN for a single cell has also been deposited in 
BioModels Database (biomodels.org) and assigned the identi¬ 
fier MODEL 1506170000. 


Representation of the solution space of constraint-based 
metabolic networks as a convex polytope 

In mass-balance setups, a set of reaction fluxes f = {fr} 
(r = 1,... ,N) describes a non-equilibrium stationary state if 
it satisfies the conditions 

N 

J2Srmfr = 0 , m = 1, . . . ,M , (29) 

r—1 

with prescribed bounds /^T < /r < /^. (Clearly, the zero 
flux vector is excluded from these considerations.) Here, 
Srm denotes the stoichiometric index of chemical species m 
{m = 1,... M) in reaction r. From a geometric viewpoint, 
( |29| l defines a convex polytope, while every point inside it 
represents a feasible flux configuration. In absence of criteria 
allowing to select specific flux vectors from the solution space 
(e.g. via the maximization of an objective function), uniform 
sampling provides important information about average fluxes 
and correlations that not only describe the viable flux patterns 
in statistical terms, but may also give relevant guidelines for 
designing experiments ll52l . When the dimension of the so¬ 
lution space is sufficiently small (typically around 10), Monte 
Carlo sampling, including straightforward rejection methods 
ES, can be applied. At higher dimensions, instead, exceeding 
computing times require the use of more advanced techniques. 


Dimensional reduction of the solution space 

We are interested in sampling the solution space of ( [29] l 
with S = {Sr-m} given by the HCCN with the prescribed 
probability measure defined in ( |26] l. Mass balance conditions 
constrain many of the 75 fluxes to be dependent on other vari¬ 
ables. To explore efficiently the space of viable flux config¬ 
urations, it is convenient to solve all such dependencies ana¬ 
lytically and then sample the (much smaller) space spanned 
by independent fluxes. To begin with, we transformed S to its 
Reduced Row Echelon Form (RREF) through Gauss-Jordan 
elimination , which can be implemented by standard pack¬ 
ages. Because the RREF of any matrix is unique, the RREF 
of the stoichiometric matrix applied to the reaction network 
uniquely represents the dependent fluxes as a function of the 
independent ones. The single-HCCN model turns out to have 
only 17 independent variables, while the two-HCCN version 
has 34. These numbers also represent the dimensions D of the 
respective solution space polytopes. 


Hit-and-run method allows to sample convex polytopes 
efficiently 

The Hit-and-Run (HR) method ||54l |55] can uniformly and 
efficiently sample a convex polytope, provided one starts from 
a point inside it. The problem of finding such a point can be 
easily solved, e.g., via Motzkin relaxation lISTl l5^ . In brief, 
HR builds a Markov chain in two steps. First, starting from 
a point inside the polytope, a random direction is extracted. 
Second, along this direction a new point is chosen uniformly 
at random inside the polytope. The scheme is then iterated 
starting from the new internal point. HR is efficient because 
once the random direction is generated the search space is re¬ 
duced to a segment. In our implementation, we generate the 
random direction using the Marsaglia-Bray method ll57l . To 
pick a point at random along the segment, instead, we first 
compute the two intersections with the polytope and then ex¬ 
tract a uniform random number inside the segment identified 
by these two points on the line. 


The Lovasz method speeds up the sampling of heterogeneous 
polytopes 

The above procedure, while fully exact, may however 
present a drawback in concrete cases. In fact, the mixing 
time of HR in convex domains scales as 0[D'^{R/r)'^] Monte 
Carlo steps, where D is the spatial dimension of the poly¬ 
tope while r and R are, respectively, the radius of the max¬ 
imum inscribed hypersphere and of the minimum circum¬ 
scribed hypersphere ifSSll . It is clear that the R/r factor can 
increase convergence times dramatically for strongly hetero¬ 
geneous polytopes. This is precisely the case for the space 
of feasible steady states in a large-scale metabolic network, 
where flux distributions can be so heterogeneous as to span 
5 orders of magnitudes lEMIl- For the network analyzed 
here, the ratio between the ellipsoid’s longest and smallest 
axes turned out to be of order 10^, leading to a condition 
number of order ~ 10®, a value that would render 

straightforward HR too expensive in terms of CPU time. To 
circumvent this ill-conditioning problem, we have resorted to 
a standard pre-processing step that identifies the ellipsoid that 
best approximates the shape of the polytope ll62ll . Extracting 
vectors uniformly on such an ellipsoid guarantees that the di¬ 
rections along the longer axis are chosen more often, thereby 
removing ill conditioning. In quantitative terms, such a pre¬ 
processing gives HR a mixing time of 0[D'^/'^] Monte Carlo 
steps, thus implying a fully polynomial scaling with D. After 
pre-processing the slowest HR mode was found to decorre¬ 
late in an affordable time of 0(10^) Monte Carlo steps. Fur¬ 
thermore, the overall pre-processing time was negligible com¬ 
pared to the time spent for the actual solution space sampling, 
while the center of the ellipsoid provides an excellent starting 
point for the algorithm. For more details of this method we 
refer the reader to ll63l . 
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A method for tuning the optimization of a given metaholic flux 
for constrained based metaholic networks 

HR is easily modified to sample points according to any 
given flux distribution P(f) by imposing that the selection 
of points along the segment takes place according to how 
the distribution P projects onto the chosen direction rather 
than uniformly. To modulate P via a linear function L{i) of 
the fluxes, it is convenient to use the Boltzmann distribution 
P(f) oc exp [/3L(f)], where /3 > 0 is an interpolation pa¬ 
rameter. By varying /3, one passes from an uniform sampling 
(corresponding to /3 = 0) to sampling flux configurations that 
maximize L(f) (corresponding to /3 ^ 1). For the HCCN 
network studied here, the ATP production /atp is maximized 
using the above method with P(f) cx exp [/3/atp]- The value 
of P for which /atp is effectively maximized can be deter¬ 
mined numerically. In Fig. S 6 we plot the ATP production as 
a function of the /3 for five different values of the maximal 
glucose available for two coupled HCCN cells. One sees that 
optimal /atp is already achieved when (3 ~ 50. 

To produce the graphs presented in Figs, [^through and 
S1 through S 6 , we applied the Lovasz method to the reduced 
row echelon forms of the HCCN networks and sampled the 


corresponding regularized polytopes by means of the Hit-and- 
run method. For curves with /3 7 ^ 0, we used the Boltzman 
weight presented above. The apply the Lovasz method and 
the uniform and optimized HR sampling, we used C-H- pro¬ 
grams that we developed and that can be freely downloaded at 
http;//chimera.romal.infn.it/SYSBIO 
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SUPPLEMENTARY TEXT 

Text SI. Human Catabolic Core Network: biosynthetic reactions 

Besides the major central carbon metabolic pathways, i.e. Glycolysis, Pentose Phosphate Pathway, Citric Acid Cycle, and 
Oxidative Phosphorylation, our model includes most reactions of Glutamate metabolism as well as lumped reactions describing 
the biosynthesis of all non-essential amino acids (i.e. the amino acids that can be synthesized by human cells) and palmitate, a 
key precursor of other fatty acids. 

Amino acid biosynthetic processes can be divided roughly in two classes. The synthesis of tyrosine and cysteine requires 
respectively the essential amino acids phenylalanine and methionine. All other amino acids can instead be synthetized from 
glutamine plus metabolites from the central carbon pathways. The processes leading to amino acid synthesis are of various 
types, from simple single transaminase reactions to a complex chain of enzymatic reactions. In the latter case, by means of mass 
balance arguments we can lump together individual reaction steps to obtain effective reactions whose net output corresponds to 
the desired amino acid(s). In all cases, our starting points were the reactions and pathways included in the reference metabolic 
network reconstruction Recon-1 [21]. In the . xml file provided as Supporting Material, we include the list of the corresponding 
Recon-1 reactions that constitute the full pathway. In the following, we briefly outline these effective reactions. 

• Glutamate synthesis: from glutamine via glutaminase (GLUN): 

gln_L -f H 2 O ^ glu_L -f NH 4 (30) 

• Alanine and aspartate synthesis: from glutamate through the transaminase reactions 

akg-I-ala_L O glu_L-I-pyr (31) 

akg -I- asp_L o glu_L -I- oaa (32) 

(akg = Qf-ketoglutarate ; pyr = pyruvate ; oaa = oxaloacetate) 

• Proline synthesis: from glutamate via glutamate 5-kinase (GLU5K), glutamate-5-semialdehyde dehydrogenase (G5SD), 
L-glutamate 5-semialdehyde dehydratase (G5SAD) and pyrroline-5-carboxylate reductase (P5CR): 

gluX -f ATP -P (2)NADPH + (2)H+ ^ proX -P ADP -P Pi + H 2 O -P (2)NADP (33) 

• Serine synthesis: from glutamate via phosphoglycerate dehydrogenase (PGCD), phosphoserine transaminase (PSERT) 
and phosphoserine phosphatase (PSP): 

3pg -P NAD -P H 2 O + gluX ^ serX -P akg -P NADH -P Pi + H+ (34) 

(3pg = 3-phosphoglyceric acid) 

• Glycine synthesis: from serinevia glycine hydroxymethyltransferase (GHMT), methylenetetrahydrofolate dehydrogenase 
(MTHFD), methenyltetrahydrofolate cyclohydrolase (MTHFC) and formyltetrahydrofolate dehydrogenase (FDH): 

serX -P H 2 O -P (2)NADP ^ gly -P CO 2 + (2)NADPH + (2)H+ (35) 

• Arginine synthesis: from glutamate and aspartate via glutamate 5-kinase (GLU5K), glutamate-5-semialdehyde dehydroge¬ 
nase (G5SD), ornithine transaminase (ORNTA), carbamoyl-phosphate synthase (GBPS), ornithine carbamoyltransferase 
(OCBT), argininosuccinate synthase (ARGSS) and argininosuccinate lyase (ARGSF): 

aspX -P gluX -P glnX -P (5)ATP -P NADPH + ( 3 )H 20 -P CO 2 ^ 

^ argX -P fum -P akg -P (5)ADP + (5)Pi -P NADP + (5)H+ (36) 

• Asparagine synthesis: from aspartate and glutamine via asparagine synthase (ASNS): 

aspX -P (2)ATP -P glnX -P ( 2 )H 20 -> (2)ADP -P asnX -P glu -P (2)H+ + (2)Pi (37) 

• Cysteine synthesis: from serine and methionine via methionine adenosyltransferase (METAT), adenosylhomocys- 
teinase (AHC), cystathionine beta-synthase (CBS), cystathionine gamma-lyase (CSE), 2-Oxobutanoate dehydrogenase, 
Propionyl-CoA carboxylase (PCCA) and methylmalonyl-CoA mutase (MMM): 

metX -P serX -P coa -P (4)ATP -P NAD + ( 4 )H 20 —> 

^ cysX -P succoa -P (4)ADP + (4)Pi -P (4)H+ + NADH -P NH 4 -P CH 3 (38) 
(coa = coenzyme A ; succoa = succinyl-coa) 
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• Tyrosine synthesis: from phenylalanine via L-Phenylalanine hydroxylase, tetrahydrobiopterin oxidoreductase, 
tetrahydrobiopterin-4a-carbinolamine dehydratase and dihydropteridine reductase: 

phe_L + NADH + O 2 + H+ ^ tyr_L + NAD + H 2 O (39) 


• Palmitic acid synthesis: 

(8)accoa + (23)ATP + (14)NADPH + (17)H20 ^ 

^ palmitate + (8)coa + (14)NADP + (23)ADP + (23)Pi + (10)H+ (40) 


(accoa = acetyl-coa) 

In addition, we included the reactions for superoxyde reduction and FADH 2 oxidation, that are necessary to correctly account 
for the redox state of the cell. For sakes of simplicity, we lumped together the reactions catalyzed by superoxide dismutase, 
gluthatione peroxidase and glutathione reductase to the single reaction 

NADPH + + 2H+ ^ NADP + 2 H 2 O . (41) 

Likewise, we unified the two reactions that oxidate FADH 2 , thereby reducing ubiquinone by means of the Electron transport 
flavoprotein-ubiquinone oxidoreductase (ETF) to 

FADH 2 + QIO ^ FAD + QIOH 2 . (42) 


Text S2. Maximizing the flux of precursors: results 

In the main text, we show that lactate overflow and shuttle arises under both biomass and ATP maximization. Here, we show 
the metabolic flux patterns obtained when in the lactate donor the production of single biomass precursors are maximized, to 
assess the robustness of the lactate overflow and shuttle. To limit the amino acid productions, we have assumed a fixed maximum 
total influx of glutamine, phenylalanine and methionine (equals to 2 mmol/(gDWh)) and a variable overall glucose supply. 

All objective functions we tested yield a lactate overflow with lactate shuttle towards the acceptor, suggesting that the crowding 
constraint is responsible for the effect. In particular, we recover high levels of lactate shuttling for palmitate optimization and 
moderate levels for the optimization of the production of amino acids. The only lactate shuttle truly related to an energetic 
imbalance is the one induced by palmitate optimization, while the others are similar in magnitude to the shuttling present in 
the absence of any optimization (i.e. for /? = 0). In Fig. S2, we display the production flux and lactate shuttling profiles as a 
function of the total glucose supply for the two-cell system when the lactate donor maximizes palmitate, proline, cysteine, and 
serine, respectively. Other biosynthetic objective functions lead to similar features. 
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SUPPLEMENTARY FIGURES 



FIG. SI. Glucose intake, ATP production, and oxidative and fermentative fluxes in a donor-acceptor symmetric system, (a, top left) 
Glucose intakes for two coupled symmetric cells as a function of the total glucose available to the donor-acceptor pair, (b, top right) ATP 
produced by the donor and the acceptor cells as a function of the total glucose available to the donor-acceptor pair, (c, bottom left)-(d, bottom 
right) Average flux through LDH (circles) and PDHm (crosses) as a function of the average glucose supplied to the two-cell system. Curves 
describe the behaviour obtained for two coupled symmetric HCCN cells with an ATP-maximizing donor (black lines, /3 = 50) or for an 
unbiased sampling of the two-cell solution space (blue lines, /3 = 0). Error bars, which represents s.e.m., are smaller than symbol sizes. 





FIG. S2. Lactate shuttling for alternative objective function maximizations. In blue we plot the flux that is maximized in the donor cell: 
palmitate (top left), proline (top right), cysteine (bottom left), and serine (bottom right) as a function of the total glucose supply, for a system 
formed by a lactate donor and a lactate acceptor. The lactate fluxes of donor and acceptor cells are depicted in green and red, respectively. A 
negative flux correspond to lactate influx. 
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FIG. S3. Pearson correlation coefficients among all fluxes of lactate donor and lactate acceptor for maximum glucose uptake Ua = 0.3. 

The fluxes of the acceptor cell are identified by the trailing letters “acc”. 





^ isiil. 



FIG. S4. Pearson correlation coefficients among all fluxes of lactate donor and lactate acceptor for maximum glucose uptake Ug = 1- 

The fluxes of the acceptor cell are identified by the trailing letters “acc”. 
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FIG. S5. Pearson correlation coefficients among all fluxes of lactate donor and lactate acceptor for maximum glucose uptake Ua = 3. 

The fluxes of the acceptor cell are identified by the trailing letters “acc”. 
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FIG. S6. By increasing /3, an HCCN increases the ATP production eventually saturating the capacity at around P ~ 50. We display 
normalized ATP production fluxes for maximal glucose supply Umax equals to 0.5 (black circels), 1 (red squares), and 5 (blu diamonds). 
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SUPPLEMENTARY TABLES 


Variable symbol 

Flux identified 

Ug 

Glucose supply 

/glyc, /hex 

Glycolysis fiux 

fox, /pDH 

Oxidative phosphorylation flux 

/ldh 

Flux through lactate dehydrogenase 

/atp 

ATP production 


TABLE S1. List of more relevant variables appearing in the mannscrlpt. 
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Metabolite short name 

Metabolite extended name 

13DPG 

3-Phospho-D-glyceroyl phosphate 

23 DPG 

2,3-Disphospho-D-glycerate 

2PG 

D-Glycerate 2-phosphate 

3PG 

3-Phospho-D-glycerate 

6 PGC 

6 -Phospho-D-gluconate 

6 PGL 

6 -phospho-D-glucono-1,5-lactone 

ACCoA 

Acetyl-CoA 

ADP 

Adenosine diphosphate 

AKG 

2-Oxoglutarate 

ATP 

Adenosine triphosphate 

CIT 

Citrate 

CO2 

Carbon dioxyde 

CoA 

Coenzyme A 

DHAP 

Dihydroxyacetone phosphate 

E4P 

D-Erythrose 4-phosphate 

F 6 P 

D-Fructose 6 -phosphate 

FAD 

Flavin adenine dinucleotide oxidized 

FADH 2 

Flavin adenine dinucleotide reduced 

FDP 

D-Fructose 1,6-bisphosphate 

FICYTC 

Ferricytochrome C 

FOCYTC 

Ferrocytochrome C 

FUM 

Fumarate 

G3P 

Glyceraldehyde 3-phosphate 

G 6 P 

D-Glucose 6 -phosphate 

GLC 

D-Glucose 

H2O 

Water molecule 

H[M] 

Hydrogen ion as an electromotive force (mitochondrial) 

H 

Hydrogen ion in cytoplasm 

ICIT 

Isocitrate 

LAC-L 

L-Lactate 

MAL-L 

L-Malate 

NAD 

Nicotinamide adenine dinucleotide 

NADH 

Nicotinamide adenine dinucleotide - reduced 

NADP 

Nicotinamide adenine dinucleotide phosphate 

NADPH 

Nicotinamide adenine dinucleotide phosphate - reduced 

02 

Molecular oxygen 

O 2 S 

Superoxide anion 

OAA 

Oxaloacetate 

PFP 

Phosphoenolpyruvate 

Pi 

Phosphate 

PYR 

Pyruvate 

QIO 

Ubiquinone-10 

QIOH 2 

Ubiquinol-10 

R5P 

alpha-D-Ribose 5-phosphate 

RU5P-D 

D-Ribulose 5-phosphate 

S7P 

Sedoheptulose 7-phosphate 

SUCC 

Succinate 

SUCCoA 

Succinyl-CoA 

XU5P-D 

D-Xylulose 5-phosphate 

giiii 

Glutamine 

glU-L 

Glutamate 

NH 4 

Ammonia 

CH 3 

Methil group 

asp_L 

Aspartate 

ala_L 

Alanine 

asn_L 

Asparagine 

pro_L 

Proline 

ser_L 

Serine 

giy 

Glycine 

argX 

Arginine 

cysX 

Cysteine 

metX 

Methionine 

tyr_L 

Tyrosine 

phe_L 

Phenylalanine 

hdca 

Palmitic acid 


TABLE S2. List of metabolites appearing in a single HCCN model. 
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Enzyme 

Reaction 

ACONT 

CIT ^ ICIT 

ACYP 

13 DPG + H 2 O ^ 3 PG + H + Pi 

AKGDm 

AKG + CoA + NAD ^ COj + NADH + SUCCoA 

ATPS4 m 

ADP + 4H + Pi ^ ATP + 3H[M] + HjO 

CSm 

ACCoA + H 2 O + OAA ^ CIT + CoA + H[M] 

CYOOm3 

4focytC + 7.92H[M] + 02 ^ 4ficytC + 4H + I. 96 H 2 O + O.O 2 O 2 

CYOR_ulOm 

2 ficytC + 2 H[M] + Q 10 H 2 ^ 2 focytC + 4 H + 

DPGM 

13DPG^ 23DPG + H 

DPGase 

23 DPG + H 2 O ^ 3 PG + Pi 

ENO 

2PG ^ H 2 O + PEP 

EBA 

EDP ^ DHAP + G3 P 

EUM 

EUM + H 2 O ^ MAL-L 

G6PDH2r 

G 6 P + NADP ^ 6 PGL + H + NADPH 

GAPD 

G3 P + NAD + Pi ^ 13 DPG + H + NADH 

GND 

6 PGC + NADP CO 2 + NADPH + RU5 P-D 

HEXl 

ATP + GLC ^ ADP + G 6 P + H 

ICDHxm 

ICIT + NAD ^ AKG + CO 2 + NADH 

ICDHy 

ICIT + NADP ^ AKG + CO 2 + NADPH 

LDH 

LAC-L + NAD ^ H + NADH + PYR 

MDH 

MAL-L + NAD ^ H + NADH + OAA 

NADH2_ul0m 

5H + NADH + Q10 4H + NAD + QIOH 2 

PDHm 

CoA + NAD + PYR ^ ACCoA + CO 2 + NADH 

PEK 

ATP + F 6 P ^ ADP + EDP + H 

PGI 

G 6 P ^ F 6 P 

PGK 

3 PG + ATP ^ 13 DPG + ADP 

PGL 

6 PGL + H 2 O 6 PGC + H 

PGM 

2PG ^ 3PG 

PYK 

ADP + H + PEP ^ ATP + PYR 

RPE 

RU5 P-D ^ XUS P-D 

RPI 

R5 P ^ RU5 P-D 

SUCDlm 

FAD + SUCC ^ FADH 2 + EUM 

SUCOASm 

ATP + CoA + SUCC ^ ADP + Pi + SUCCoA 

TALA 

G3P + S7P^E4P + F 6 P 

TKTl 

R5 P + XUS P-D ^ G3 P + S7 P 

TKT2 

E4P + XUSP-D^ F 6 P + G3P 

TPl 

DHAP ^ G3 P 

02S reduction 

NADPH + O 2 + 2 H -> NADP + 2 H 2 O 

PAD regeneration 

QIO + FADH 2 ^ QIOH 2 + FAD 

ATP consumption 

ATP + H 2 O ^ ADP + Pi + H 


TABLE S3. Internal reactions of an HCCN with the corresponding enzyme that catalyzes it: main catabolic pathways. 
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Enzyme 

Reaction 

ALATAX 

AKG -f alaX ^ gluX -f PYR 

ASPTA 

AKG + aspX ^ gluX + OAA 

GLUDxm 

gluX + H2O + NAD ^ AKG + H + NADH -f NH4 

GLUDym 

gluX + H2O + NADP ^ AKG + H -f NADPH + NH4 

GLUN 

glnX H2O gluX -b NH4 

Tyrosine production 

pheX + NADH + O2 -b h —tyrX -b NAD -b H2O 

Asparagine production 

aspX -b (2)ATP + glnX -b (2)H20 —> (jfADP + asnX + gluX -b (2)H + (2)Pi 

Proline production 

gluX -b ATP -b (2)NADPH -b 2 H^ proX -b ADP -b Pi -b HjO -b (jiNADP 

Serine production 

3 PG + NAD -b H2O + gluX ^ serX -b AKG + NADH + Pi + H 

Glycine production 

serX -b H2O -b (,)NADP gly -b CO2 + (2)NADPH -b (2)H 

Arginine production 

aspX + gluX + glnX -b (jiATP + NADPH -b 3 H2O + CO2 ^ 
argX + EUM + AKG + (5)ADP -b (jfPi + 5 H + NADP 

Cysteine production 

metX -b serX + GOA -b (4)ATP + NAD -b CO2 -b (4)H20 
cysX -b SUCCOA + (4)ADP + (4)Pi -b 4 H -b NADH -b NH4 -b CO2 -b CH3 

Palmitic acid production 

(gfACCOA -b (23)ATP + 8 NADH -b (s)NADPH -b 17 H2O ^ 
hdca -b (sfCOA -b (sfNAD -b (sfNADP -b (23)ADP + (23)Pi + (,o)H 


TABLE S4. Internal reactions of an HCCN with the corresponding enzyme that catalyzes it: glutamate metabolism and anabolic 
effective reactions. 


Description 

Lower bound 

Upper bound 

Reaction 

Flux of carbon dioxide 

-1000 

1000 

C 02 -S—out 

Flux of oxygen 

-1000 

1000 

O 2 <—> out 

Flux of water 

-1000 

1000 

H 2 O ^out 

Flux of ion hydrogen 

-1000 

1000 

H <—> out 

Glucose intake 

0 

Ua 

in ^— GLC 

Lactate flux 

0 

-1000 

1000 

1000 

LAC-L —out (donor) 
LAC-L <—> out (acceptor) 

Glutamine intake 

0 

1 

in <— glnX 

Phenylalanine intake 

0 

1000 

in <— pheX 

Methionine intake 

0 

1000 

in <— metX 

Ammonia flux 

-1000 

1000 

NH 4 out 

Methyl group flux 

-1000 

1000 

CH 3 out 


TABLE S5. Exchange reactions of an HCCN with the corresponding bounds. In this model, we assume that carbon dioxide, water, oxygen, 
and proton can freely diffuse in and out of the cell. The maximum glucose intake for each cell cannot exceed the total glucose supply. Lactate 
flux distinguishes donor and acceptor cells as the donor can only secrete lactate, while the acceptor can also intake lactate. In Table [S^ the 
constraints on glucose and lactate fluxes for the donor-acceptor couple are presented. 
































25 


Enzyme 

Lower bound 

Upper bound 

Enzyme extended name 

ACONT 

-1000 

1000 

Aconitase 

ACYP 

0 

1000 

Acylphosphatase 

AKGDm 

0 

1000 

2-Oxoglutarate dehydrogenase 

ATPS4m 

0 

1000 

ATP synthase 

CSm 

0 

1000 

Citrate synthase 

CYOOm3 

0 

1000 

Cytochrome C oxidase, mitochondrial Complex IV 

CYOR_ulOm 

0 

1000 

Ubiquinol-6 cytochrome C reductase, Complex III 

DPGM 

-1000 

1000 

Diphosphoglyceromutase 

DPGase 

0 

1000 

Diphosphoglycerate phosphatase 

ENO 

-1000 

1000 

Enolase 

FBA 

-1000 

1000 

Fructose-bisphosphate aldolase 

FUM 

-1000 

1000 

Fumarase 

G6PDH2r 

-1000 

1000 

Glucose 6-phosphate dehydrogenase 

GAPD 

-1000 

1000 

Glyceraldehyde-3-phosphate dehydrogenase 

GND 

0 

1000 

Phosphogluconate dehydrogenase 

HEXl 

0 

1000 

Hexokinase 

ICDHxm 

0 

1000 

Isocitrate dehydrogenase 

ICDHy 

0 

1000 

Isocitrate dehydrogenase 

LDH 

-1000 

1000 

Lactate dehydrogenase 

MDH 

-1000 

1000 

Malate dehydrogenase 

NADH2_ul0m 

0 

1000 

NADH dehydrogenase, mitochondrial 

PDHm 

0 

1000 

Pyruvate dehydrogenase 

PFK 

0 

1000 

Phosphofructokinase 

PGI 

-1000 

1000 

Glucose-6-phosphate isomerase 

PGK 

-1000 

0 

Phosphoglycerate kinase 

PGL 

0 

1000 

6 -phosphogluconolactonase 

PGM 

-1000 

1000 

Phosphoglycerate mutase 

PYK 

0 

1000 

Pyruvate kinase 

RPE 

-1000 

1000 

Ribulose-5-phosphate 3-epimerase 

RPI 

-1000 

1000 

Ribose-5-phosphate isomerase 

SUCDlm 

-1000 

1000 

Succinate dehydrogenase 

SUCOASm 

-1000 

1000 

Succinate-CoA ligase 

TALA 

-1000 

1000 

Transaldolase 

TKTl 

-1000 

1000 

Transketolase 

TKT2 

-1000 

1000 

transketolase 

TPI 

-1000 

1000 

Triose-phosphate isomerase 

Lumped reaction 

0 

1000 

Reduction of superoxyde anion 

Lumped reaction 

0 

1000 

FAD regeneration 

Effective reaction 

T*min 

/ATP 

1000 

ATP hydrolysis 


TABLE S6. Internal reactions of an HCCN with the corresponding hounds: catabolic pathways. 


Enzyme 

Lower bound 

Upper bound 

Enzyme extended name 

ALATA_L 

-1000 

1000 

alanine transaminase 

ASPTA 

-1000 

1000 

aspartate transaminase 

GLUDxm 

-1000 

1000 

glutamate dehydrogenase (NAD) 

GLUDym 

0 

1000 

glutamate dehydrogenase (NADP) 

GLUN 

0 

1000 

glutaminase 

Lumped reaction 

0 

1000 

Tyrosine production 

Lumped reaction 

0 

1000 

Asparagine production 

Lumped reaction 

0 

1000 

Proline production 

Lumped reaction 

0 

1000 

Serine production 

Lumped reaction 

0 

1000 

Glycine production 

Lumped reaction 

0 

1000 

Arginine production 

Lumped reaction 

0 

1000 

Cysteine production 

Lumped reaction 

0 

1000 

Palmitic acid production 


TABLE S7. Internal reactions of an HCCN with the corresponding hounds: glutamate metabolism and anabolic effective reactions. 
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Variable constrained 

Lower bound 

Upper bound 

Variable expressed by reactions 

Donor resource 

0 

4? ATP 

fflglyc/HEX.don + fflox/pDH.don + aLDn/LDH.don 

Acceptor resource 

0 

4? ATP 

fflglyc/HEX.acc + flox/pDH.acc + OLDH ]/pDH.acc | 

Total glucose intake 

0 

Ug 

^G,don ^G,acc 

Total lactate flux 

0 

1000 

^LAC,don H” t^LAC,acc 


TABLE S8. Constraints on maximum resources available for ATP production and on total glucose and lactate fluxes. The first column 
contains a description of the variable that is constrained, while the last column the variable expressed as a function of the elementary fluxes of 
the metabolic network. The second and the third column contain the minimum and the maximum value that the variable can take, respectively. 
<1?ATP is set to 0.4 and states that at most 40% of the cellular resources can be devoted to ATP production. The sum of the glucose intaken by 
donor and acceptor cannot exceed the total glucose supply Ug- The constraint on total lactate flux establishes that there is no external source 
of lactate in the system and that the donor-acceptor couple can only produce lactate. 


Metabolite 

Biomass coefficient 

H 2 O 

-20.651 

ATP 

-20.651 

ADP 

20.651 

H 

20.651 

Pi 

20.971 

glu-L 

-0.38587 

asp_L 

-0.35261 

asn_L 

-0.27942 

ala_L 

-0.50563 

cysX 

-0.046571 

gln-L 

-0.326 

giy 

-0.53889 

ser_L 

-0.39253 

arg_L 

-0.35926 

met_L 

-0.15302 

tyr_L 

-0.15967 

phe_L 

-0.25947 

pro_L 

-0.41248 

hdca 

-0.112 

R5P 

-0.045 

G6P 

-0.275 


TABLE S9. Biomass objective function coefficients adapted from [47]. 




